function pcut = addBilinearVariableCuts(p)

if isempty(p.bilinears)
    pcut = p;
else
    pcut = emptyNumericalModel;
    
    z = p.bilinears(:,1);
    x = p.bilinears(:,2);
    y = p.bilinears(:,3);
    
    nn = length(p.c);
    
    still_uncertain = find(abs(p.lb(z)-p.ub(z))>1e-8);
    if ~isempty(still_uncertain)
        
        x_lb = p.lb(x);
        x_ub = p.ub(x);
        y_lb = p.lb(y);
        y_ub = p.ub(y);
        m = length(x);
        one = ones(m,1);
        general_vals =[x_lb.*y_lb one -y_lb -x_lb,x_ub.*y_ub one -y_ub -x_ub,-x_ub.*y_lb -one y_lb x_ub,-x_lb.*y_ub -one y_ub x_lb]';
        general_cols = [one z+1 x+1 y+1 one z+1  x+1 y+1 one z+1 x+1 y+1 one z+1 x+1 y+1]';
        general_row = [1;1;1;1;2;2;2;2;3;3;3;3;4;4;4;4];
        
        quadratic_row = [1;1;1;2;2 ;2;3; 3; 3;4;4;4;5;5;5;6;6;6];
        quadratic_cols =  [one  z+1 x+1 ,one z+1 x+1 ,one  z+1 x+1, one z+1 x+1 ,one z+1 x+1, one z+1 x+1]';
        x1 = (3*x_ub+x_lb)/4;
        x2 = (x_ub+3*x_lb)/4;
        x3 = (x_ub+x_lb)/2;
        quadratic_vals = [-x_ub.*x_lb -one x_lb+x_ub x_lb.*y_lb one -y_lb-x_lb x_ub.*y_ub one -y_ub-x_ub  x1.*x1 one -x1-x1 x2.*x2 one -x2-x2 x3.*x3 one -x3-x3]';
        m = 1+length(p.c);
        rows = [];
        cols = [];
        vals = [];
        nrow = 0;
        
        for i =still_uncertain(:)'
            x = p.bilinears(i,2);
            y = p.bilinears(i,3);
            if x~=y
                rows = [rows;general_row+nrow];
                vals = [vals;general_vals(:,i)];
                cols = [cols;general_cols(:,i)];
                nrow = nrow + 4;
            else
                col = quadratic_cols(:,i);
                val = quadratic_vals(:,i);
                
                rows = [rows;quadratic_row+nrow];
                vals = [vals;val];
                cols = [cols;col];
                nrow = nrow + max(quadratic_row);
            end
        end
        
        F_temp = sparse(rows,cols,vals,nrow,m);
        keep = find(~isinf(F_temp(:,1)) & ~isnan(F_temp(:,1)));
        F_temp = F_temp(keep,:);
        pcut.F_struc = [F_temp;pcut.F_struc];
        pcut.K.l = pcut.K.l+size(F_temp,1);
    end
    
    pcut = mergeNumericalModels(p,pcut);
end



